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We present our numerical results of two-dimensional hydrodynamic (HD) simulations and magnetohydrodynamic 
(MHD) simulations of the collapse of rotating massive stars in light of the collapsar model of gamma-ray bursts 
(GRBs). Pushed by recent evolution calculations of GRB progenitors, we focus on lower angular momentum of 
the central core than the ones taken mostly in previous studies. By performing special relativistic simulations 
including both realistic equation of state and neutrino cooling, we follow a long-term evolution of the slowly 
rotating collapsars up to ~ 10 s, accompanied by the formation of jets and accretion disks. We find such 
outfows can be launched both by MHD process and neutrino process. We invest igeate the properties of these 
jets whether it can become GRBs or remains primary weak outflow. 



1. Introduction 

Gamma-ray bursts (GRBs) are one of the most en- 
ergetic phenomena in the universe. Pushed by recent 
observations, the so-called collapsar has received quite 
some interest for the central engines of the long GRBs 
d, QjJ EH • Here we focus on the outflow formation in 
collapsars. There are two promising senario to launch 
outflow from collapsar, MHD process and/or neutrino 
process, which has also been extensively investigated 
thus far (e.g., @, i, g 1, M El El ) • 

A general outcome of the MHD studies among them 
is that if the central progenitor cores have significant 
angular momentum (~ 10 17 cm 2 /s) with strong mag- 
netic fields (~ 10 11 G), magneto-driven jets can be 
launched strong enough to expel the matter along the 
rotational axis within several seconds after the onset 
of collapse. It should be noted that the angular mo- 
mentum of those GRB progenitors (~ 10 16 cm 2 /s), 
albeit not a final answer due to much uncertainty, is 
relatively smaller than those assumed in most of the 
previous collapsar simulations. These situations moti- 
vate us to focus on slower rotation of the central core 
in the context of collapsar models Q . As for the initial 
magnetic fields, we choose to explore relatively smaller 
fields (< 10 10 G), which has been less investigated so 
far. Paying particular attention to the smaller angu- 
lar momentum, it takes much longer time to amplify 
the magnetic fields large enough to launch the MHD 
jets than previously estimated. By performing long- 
term evolution over ~ 10 s, we aim to clarify how the 
properties and the mechanism of the MHD jets could 
change with the initial angular momentum and the 
initial magnetic fields. 

In terms of the neutrino-driven outflows, it still re- 
mains preliminary ones, such as evaluation in the post 



processing manner [2j] . We have developed code to cal- 
culate neutrino heating via neutrino-antineutrino pair 
annihilation with ray-tracing method [6|. We show 
some result of the test calculation and application to 
the collapsar model. 



2. Magneto-Driven Jet from Collapsars 
2.1. Numerical Methods and Initial Model 

The results presented here are calculated by the 
MHD code in special relativity developed by [5|, [l5| , 
including both realistic equation of state and neu- 
trino cooling. As for the initial profiles of the col- 
lapsing star, we employ the spherical data set of den- 
sity, temperature, internal energy, and electron frac- 
tion in model 350C in [17|. We add cylindrical ro- 
tation and dipole magnetic field profiles with vanish- 
ing toroidal magnetic field in a parametric manner. 
We compute 15 models changing the initial rotation 
and the strength of magnetic fields (See [5[ for de- 
tails). Model name is given as 'BXJY' representing 
the model with B = 10 X G and j — Yj\ 80 , where j is 
the specific angular momentum and ji so is the specific 
angular momentum for the last stable angular orbit. 



2.2. Result 

Computing 15 models in a longer time stretch than 
ever among previous collapsar models, we observe a 
wide variety of the dynamics changing drastically with 
time. Here we summarize main results of this study. 

To extract the general features furthermore among 
the models, we focus on the properties of MHD out- 
flows and neutrino luminosities. It is noted that both 
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Model 


BIO 


B9 


B8 




O (TYPE II) 


O (TYPE I) 


X 






1.1 x 10 52 erg/s 


1.6 x 10 52 erg/s 


1.8 x 10 52 


erg/s 


J1.5 


3.4 s 


9.2 s 


4.3 s 






O (TYPE II) 


O (TYPE I) 


X 






4.5 x 10 51 erg/s 


5.1 x 10 51 erg/s 


8.5 x 10 51 


erg/s 


J2.0 


5.8 s 


5.3 s 


6.0 s 






O (TYPE I) 


X 


X 






1.6 x 10 50 erg/s 


1.4 x 10 50 erg/s 


1.7 x 10 50 


erg/s 


J2.5 


7.7 s 


12 s 


10 s 






O (TYPE I) 


X 


X 






9.0 x 10 49 erg/s 


2.5 x 10 49 erg/s 


4.5 x 10 49 


erg/s 


J3.0 


9.3 s 


14 s 


12 s 





-3e+07 



-1.5e+07 



1.5e+07 



3e+07 



Table I Properties of MHD jets and neutrino 
luminosities. Contents of each cell is, whether the MHD 
jets are formed (yes or no indicated by O or x ) with the 
different formation mechanisms indicated by TYPE I or 
TYPE II (top), the neutrino luminosity (middle) 
estimated at the epoch (bottom) when the accretion 
disks become almost stationary. 
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Figure 1: Logarithmic contour of differential rotation 
(:XdQ/dX) (left) and amplification rate of toroidal 
magnetic field (-.dB^/dt) (right). 



of them is helpful to understand the energy sources 
for powering the GRBs namely via magnetic and/or 
neutrino-heating mechanisms. 

Top column of each cell of Table [J indicates whether 
the MHD jets are formed (yes or no indicated by O 
or x). TYPE I or II indicates the difference of the 
formation process of the MHD jets (see below). The 
quantities of the middle column show the neutrino 
luminosity (sum of all the neutrino species v ei v ei 
and vx) estimated at the epoch when the accretion 
disks becomes stationary (bottom) (e.g., typically ~ 
4 sec). We find that the neutrino luminosities be- 
come higher for slower rotation models. This is be- 
cause the accretion disks can attain higher tempera- 
tures due to the gravitational compression. It is in- 
teresting to note that the luminosities tend to become 
smaller for strongly magnetized models with relatively 
smaller angular momentum (J < J2.5). This is mainly 
because the gravitational compression is hindered by 
the magnetic forces confined in the disks. 

2.2.1 . amplification of magnetic field and formation of 
magnetic wind 

Since the initial models investigated here are as- 
sumed to have only the poloidal fields, the key in- 
gredients for amplifying the toroidal fields are the 
compression of the poloidal fields and the efficient 
wrapping of them via differential rotations. In ad- 
dition, MRI should also play an important role, 
whose wavelength of the fastest growing mode is given 

^ A ~ 5 (^)(^)(l^!) 1/2 km. Here, 
putting the typical physical values of the disk, our nu- 
merical grids are insufficient to capture MRI at earlier 
phase when the magnetic field is weak, but can han- 



dle it in the later phase when the magnetic field gets 
stronger. In this sense, the discussion below should 
give a lower bound for the field amplification. 

Figure [T] shows the differential rotation and am- 
plification rates of toroidal magnetic field of model 
B9J1.5 at 7.99 s, just before the launch of the MHD 
outflows from the accretion disk. The density takes 
its maximum value at around 120 km in the equa- 
torial plane (left panel), in which the poloidal fields 
are strong because the higher compression is achieved 
(right panel). The white solid line in the right panel 
indicates the surface positions of the accretion disk. 
So the amplifications of the poloidal fields occur most 
efficiently inside the accretion disk. It is noted here 
that the disk is gravitationally stable because the adi- 
abatic index (7) inside the disk becomes greater than 
4/3 due to the contribution of the non-relativistic nu- 
cleon (7 = 5/3) photodissociated from the iron nuclei. 
Figure [1] shows that the amplification rates of toroidal 
magnetic field are highest also inside the disk (right), 
because the degree of the differential rotation is large 
there (left). In previous collapsar simulations assum- 
ing much larger angular momentum initially, it seems 
to be widely agreed that the differential rotation is 
a primary agent to amplify the toroidal fields. On 
the other hand, our results show that for long-term 
evolution of relatively slow rotation models, the am- 
plification of the poloidal fields by compression is pre- 
conditioned for the amplification of the toroidal fields. 

2.2.2. two types of megneto-driven outflow 

Figure [2] shows the evolution of the MHD outflow 
for model B10 J 1.5, from its initiation near from the 
inner edge of the accretion disk (top left), propaga- 
tion along the polar axis (top right), till they come 
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Figure 2: Evolution of MHD jets launched from the 
accretion disk (Type II jets). In each panel, logarithmic 
contour of entropy (left side) and the inverse of plasma 
beta (e.g., f3~l g = B 2 /8tt/p) (right side) are shown at 
1.43 s (top left), 1.87 s (top right), and 2.11 s (bottom) 
for model B10J1.5. Note the difference of the length 
scales among panels. 




Figure 3: Three dimensional plot of density with the 
magnetic field lines (silver line) for BIO J 1.5 model near 
at the moment of the shock break-out (2.11 s). Color 
contour on the two-dimensional slice represents the 
logarithmic density. 



out of the iron core (bottom). Note the difference of 
the length scales in each panel. Among the computed 
models, this model has the strongest initial magnetic 
fields with smallest angular momentum (e.g., Table 
0}. The toroidal fields can be much stronger than 
for other models, by the enhanced compression of the 
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Figure 4: Same as Figure [2] but for model B9J1.5 at 8.115 
s (top left), 8.148 s (top right), and 8.154 s (bottom), 
showing the moment of the formations of jets in type I. 



poloidal fields inside the disk. As a result, the MHD 
outflow is so strong that they do not stall, once they 
are launched. In fact, the outflow is shown to be kept 
magnetically-dominated (inverse of the plasma beta 
greater than 1, right-hand side in Figure [2]) till the 
shock break-out. As the outflow propagates rather 
further from the center (top right), the outflow begins 
to be collimated due to the magnetic hoop stresses, 
and keep their shape till the shock-breakout (Figure 

S). 

We move to discuss the jets in type I, by taking 
model B9J1.5 as an example. In this case, the mag- 
netic outflow from the accretion disk is not as strong 
as model B10J1.5, and stalls at first in the iron core 
(see butterfly-like regions colored by red in the top 
left panel (right-hand side) in Figured]). In the top 
right panel, very narrow regions near along the ro- 
tational axis are seen to be produced in which the 
magnetic pressure dominate over the matter pressure 
(colored by red in the right-hand side). Such regions 
are formed by turbulent inflows of the accreting ma- 
terial from the equator, crossing the butterfly- like re- 
gions outside the disk, to the polar regions. Such flow- 
in materials obtain sufficient magnetic amplifications 
when they approach to the rotational axis where the 
differential rotation is stronger, leading to the forma- 
tions of the MHD outflows along the rotational axis 
(bottom). 

Now we proceed to look more in detail to the prop- 
erties of the MHD jets. The jet of model B10J1.5 has 
the largest explosion energy with largest baryon load- 
ing. This is because the jet is type II as mentioned. 
Since the jet is launched rather earlier (~ 1.9s) than 
for type I, there is much material near the rotational 
axis, which makes the baryon load of jets larger for 
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Figure 5: An example of the ray-tracing (white line) of 
neutrinos for estimating the neutrino pair annihilation 
towards a given point outside the accretion disk 
(indicated by red). The central sphere represents the 
black hole (BH). 



the model. For type I jets, no systematic dependence 
of the initial angular momentum on the masses and 
the energies is found. We think that this is because 
the formation of the type I jets occurs by turbulent 
inflows as already mentioned in the previous section. 
The similarity between types I and II is that the jets 
are at most subrelativistic (0.07c for model B10J1.5) 
with the explosion energy less than 10 49 erg. While the 
ordinary GRBs require the highly relativistic ejecta, 
we speculate that thesemildly relativistic ejecta may 
be favorable for X-ray flashes (XRF). 



3. Neutrino-Driven Jet from Collapsars 
3.1. Numerical Method 

3.1 .1 . neutrino pair annihilation in collapsars; ray-tracing 
method in special relativity 

We develop a numerical scheme and code for esti- 
mating the energy and momentum transfer via neu- 
trino pair annihilation (y + v — » e~ + e + ), bearing in 
mind the application to the collapsar models of GRBs 
@. To calculate the neutrino flux illuminated from 
the accretion disk, we perform a ray-tracing calcula- 
tion in the framework of special relativity. 

Figure [5] presents our calculation concept for the 
neutrino pair annihilation in collapsars (FigurefSJ). For 
estimating the neutrino pair annihilation at a given 
target, we trace each neutrino trajectory (white lines) 
backwards till it hits to the surface of the accretion 
disk (colored by red) . We checked the accuracy of our 
code for various situations, and confirmed its validity. 



3.1 .2. two step simulation 



In this section, we show results of the application of 
newly developed code for neutrino heating Q to col- 
lapsars. First of all, we note our basic strategy. Our 
numerical procedure is divided into two steps, STEP I 
and STEP II, depending on the inclusion of neutrino 
heating, since it is numerically expensive. In STEP 
I, we follow the evolution of collapsar from the onset 
of core-collapse without neutrino heating and evalu- 
ate the effect of neutrino heating in post-processing 
step. If the neutrino heating meets the condition for 
outflow formation, we go into STEP II in which we 
restart our simulation including neutrino heating from 
the time when the heating becomes important for dy- 
namics. The criterion of this switching is whether 
Tdyn is greater than Th ea t or not, where Td yn and 
Theat is the dynamical timescale of materials and the 
timescale for materials to escape from gravitational 
field by neutrino heating, respectively. Therefore, the 
condition Td yn /Th e at > 1 means that materials can 
escape from gravitational field by neutrino energy de- 
position. Here, the dynamical timescale is given by 
Tdyn = v^Stt /16Gp, where G is the gravitational con- 
stant and p is average density derived from spherical 
mass coordinate M{r) as p = 3M(r)/47rr 3 . Heating 
time scale is given by r heat = p|$tot|/g + , where p, q + , 
and (frtot is the density, the energy deposition rate by 
neutrino heating, and the total gravitational poten- 
tial. 

In STEP I, we start our simulation from the on- 
set of core-collapse of the rotating massive star. The 
time-evolution is followed by special relativistic hydro- 
dynamic code developed by [TBI . Since the massive 
star is thought to form the BH after the onset of core- 
collapse, we simply assume the formation of the BH 
by setting absorptive inner boundary, and put a point 
mass Mbh at the center of the progenitor. 

In terms of heating, we include the energy de- 
position and the momentum transfer by neutrino- 
antineutrino pair annihilation for electron type neu- 
trino. To calculate the heating rate, we adopt the 
method in p, which is developed for the rapid calcu- 
lation of pair annihilation in collapsars. With this 
method, we evaluate the heating rate in the post- 
processing step in STEP I. Since the simulation of 
heating is still time-demanding to couple with the hy- 
drodynamics code, we calculate neutrino heating once 
in every 50 hydrodynamical steps in STEP II, which 
is enough to capture the dynamical change of the neu- 
trino flux from the accretion disk. 

We calculated 4 models with various profiles of ro- 
tation. Model name is given as 'JX' representing the 
model with angular momentum j = Yj\ so . 
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1e+50 

4 5 6 
Time [s] 

Figure 6: Time evolution of total luminosity and the 
mass accretion rate. Plots for the luminosity of JO. 6 
model (solid), JO. 8 model (dashed), J1.0 model (dotted), 
J1.2 model (small dotted), and the mass accretion rate in 
JO. 8 model (dot dashed) is shown. 
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Figure 7: Hydrodynamic configuration when the 
accretion disk is in a stationary state (here, 9.1 s after the 
onset of core-collapse). In the left panel, the logarithmic 
density (ing/cm 3 , left-half) and temperature (in K, 
right- half) are shown. The right panel shows the electron 
fraction (left-half) and the Lorentz factor (right-half). 



3.2. Result 

3.2.1 . when the neutrino heating becomes important ? 

First of all, we mention about the evolution of col- 
lapsars without neutrino heating for 11 s (STEP I). 
All of our models initially collapse into the inner ra- 
dius. When relatively rapidly rotating matter with, 
especially at the edge of Fe core, reaches near the in- 
ner radius, material crosses the equatorial plane and a 
shock is formed there (~ 0.8 s). This shock efficiently 
converts kinetic energy to thermal energy, then thick 
disk is formed in all models. At 2.0 s after the onset 
of core-collapse, Mbh becomes greater than 2 M in 
each model, which would safely verify our assumption 
of BH formation. Figure [71 shows the stationary disk 
at 9.1 s after the onset of core-collapse in JO. 8 model. 
We find that large and hot disk can be formed for 
models with angular momentum smaller than jrmiso- 

Here after we move to the effect of neutrinos emit- 
ted from the accretion disk. One of the most im- 
portant ingredients for neutrino-driven outflow would 



J0.6 


No 


5.0 x 10 50 




0.0 


J0.8 
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1.2 x 10 53 


3.1 x 10 50 


4.0 x 10 49 


J1.0 


No 


4.1 x 10 52 


2.9 x 10 49 


0.0 


J1.2 


No 


2.3 x 10 52 


8.9 x 10 48 


0.0 



Table II Properties of our models. Contents of each cell 
is, whether the neutrino-driven outflows are formed (yes 
or no), the total neutrino luminosity L u , the neutrino 
heating rate by pair annihilation Q u „, and the increasing 
rate of explosion energy E exp . 



be the total luminosity of neutrino, which will be 
potentially converted to the thermal energy. Fig- 
ure [6] shows time evolution of the neutrino luminos- 
ity. We fine that after several seconds (namely 5 
s), the neutrino luminosity can become greater than 
10 53 ergs s _1 . We also find that the criterion for out- 
flow formation Td yn /Theat > 1 is safely satisfied after 
about 9.0 s only in JO. 8 model, indicating that outflow 
would be launched. Therefore we move onto STEP II 
and restart our simulation including neutrino heating 
from 9.0 s. 

3.2.2. formation of neutrino-driven outflow 

In STEP II, we find that just after we restart with 
neutrino heating, an outflow is launched from the po- 
lar funnel in JO. 8 model, while not in other models. 
Table [TT] shows the main results of this work. In this 
table, the total luminosity L to t 5 the total energy de- 
position rate Q uP , and the luminosity of the explo- 
sion energy ^ exp of our models are shown, which are 
given by the average value from 9.1000 s till 9.2044 
s after the onset of core-collapse. As shown, E exp 
is non-zero (4.0 x 10 49 ergs s _1 ) only for JO. 8 model. 
This outflow is driven by neutrino-antineutrino pair 
annihilation, which we will discuss below. Amang 
all models, JO. 8 model has the largest luminosity 
L to t ~ 1.2 x 10 53 ergs s _1 and energy deposition rate 
Q V y ^ 3.1 x 10 50 ergs s _1 , respectively, providing the 
energy conversion efficiency of 0.26 %. Q V y is com- 
parable to the artificial energy deposition rate in pre- 
vious studies (e.g., [2, Hj]), indicating the analogeous 
outcomes in the outer layer. We note that E exp of JO. 8 
model is consistent with typical energy of 7-ray in the 
rest frame of GRBs (~ 10 51 ergs s _1 ) if this outflow 
continues for about 100 s and the conversion rate to 
kinetic energy is several 10% (e.g., H), indicating the 
connection between our model and GRBs. 



4. Summary 

We have performed a series of numerical simula- 
tion of collapsars from core-collapsing phase of mas- 
sive star. We focus on indevestigating whether outflow 
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Figure 8: Relativistic outflow production in JO. 8 model. In each panel, logarithmic contour of the entropy (left) and 
the ratio of the heating timescale to the dynamical timescale Tdyn/Theat (right) are drawn. Each panel corresponds to 
9.1000 s (left), 9.1025 s (center), and 9.1079 s (right) after the onset of core- collapse, respectively. The white solid line 
denotes the area where the density is equal to 10 n g/cm 3 , showing the surface of the accretion disk. 



can be launched by MHD and/or neutrino process. 
In terms of MHD process, we find that initial weak 
poloidal magnetic field is much amplified inside the 
disk by compression. Such strong magnetic field can 
amplify the toroidal magnetic field. We find that after 
several seconds, magnetic pressure dominates over the 
thermal pressure, and trigers the formation of outflow. 
We also find that there are two types of MHD outflows 
in long-term evolution from core-collapse. Although 
these outflows are weak (« 10 48 ergs s _1 ), it can be 
candidate for the source of XRF. 

As for the neutrino process, we developed numerical 
code for pair neutrino annihilation, and applied it to 
collapsars. We find that the condition for the outflow 
production by neutrino-antineutrino pair annihilation 
is fulfilled inside 100 km on the axis at 9 s after core- 
collapse. Including the energy deposition by neutrino, 
we find that neutrino-driven outflow is formed along 
the rotational axis. In addition, the outflow becomes 
relativistic. Such relativistic outflow can keep colli- 
mated till it penetrate the outer layer. These features 
suggests the connection between neutrino-driven out- 
flows and GRBs. 
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